# Split data based on setting vs. no setting
no_setting <- mtb_data %>% 
  filter(is.na(setting))
setting <- mtb_data %>% 
  filter(!is.na(setting))



setting <- cbind(setting$model, setting$label, select_if(setting, is.numeric))
setting$model <- setting$`setting$model`
setting <- setting %>% select(-`setting$model`)
setting$label <- setting$`setting$label`
setting <- setting %>% select(-`setting$label`)

mean_by_setting <- aggregate(x=select(setting, -c(model, label)),
                             by=list(setting$model, setting$label),
                             FUN=mean)
mean_by_setting$model <- mean_by_setting$Group.1
mean_by_setting$label <- mean_by_setting$Group.2
mean_by_setting <- mean_by_setting %>% select(-c(Group.1, Group.2))

no_setting <- cbind(no_setting$model, no_setting$label, select_if(no_setting, is.numeric))
no_setting$model <- no_setting$`no_setting$model`
no_setting <- no_setting %>% select(-`no_setting$model`)
no_setting$label <- no_setting$`no_setting$label`
no_setting <- no_setting %>% select(-`no_setting$label`)

new_mtb_data <- data.frame(rbind(mean_by_setting, no_setting))

rownames(new_mtb_data) <- new_mtb_data$model

mtb_correlation <- mtb_data %>% 
  # Get rid of price for now
  select(-price) %>% 
  # Select our variables of interest
  select_if(is.numeric) %>% 
  # Remove rows with NAs in them
  # drop_na() %>% 
  # Build our correlation matrix, such that missing values are handled by casewise deletion
  cor(use = 'complete.obs') 

# Convert our results into a tibble for easier manipulation
mtb_correlation_df <- mtb_correlation %>% 
  as_tibble() %>% 
  mutate(variable = colnames(mtb_correlation)) %>% 
  relocate(variable, everything())

# Build our correlation plot, using the angular order of the eigenvectors
corrplot(mtb_correlation,
         diag = F,
         col = COL2('PRGn'),
         tl.col = 'slateblue4',
         type = 'lower',
         method = 'color',
         order = 'AOE',
         title = 'Mountain Bike Feature Correlation'
         )

mtb_correlation_df %>% 
  pivot_longer(-variable, 
               names_to = 'correlated_variable', 
               values_to = 'correlation') %>% 
  filter(variable != correlated_variable) %>% 
  filter((correlation > .95) | (correlation < -.95)) %>% 
  pander()

# Impute missing values with column mean (not really best practice, but good enough)
for (c in 1:ncol(new_mtb_data)){
  if (is.numeric(unlist(new_mtb_data[,c]))){
    # print(colnames(new_mtb_data)[c])
    new_mtb_data[is.na(new_mtb_data[,c]), c] <- mean(unlist(new_mtb_data[,c]), na.rm=TRUE)  
  }
}

# TODO get average bikes by setting  to work above
mtb_no_null <- new_mtb_data %>% 
                select(-price) %>%
                select_if(is.numeric) %>% 
                bind_cols(label = new_mtb_data$label) %>%
                drop_na()


# head(mtb_no_null)


mtb_pca <- prcomp(mtb_no_null %>% select(-label),
                  center = TRUE,
                  scale. = TRUE)

# Put our summary results into a dataframe
mtb_pca_df <- tibble(variable = c('Standard Deviation', 'Proportion of Variance', 'Cumulative Proportion')) %>% 
  cbind(summary(mtb_pca)$importance)


# head(mtb_pca_df)

mtb_pca_df %>% 
  # Only display the first 6 columns
  select(c(variable:PC5)) %>% 
  pander()

p_load(devtools,
       ggbiplot)

ggbiplot(mtb_pca,
              obs.scale = 1,
              var.scale = 1,
              groups = mtb_no_null$label,
              ellipse = TRUE,
              circle = FALSE,
              ellipse.prob = .5) + 
  theme(legend.direction = 'horizontal',
               legend.position = 'top')

# jpeg('../Images/pca.jpg')

# How many clusters are necessary? 4

mtb_numeric <- mtb_no_null %>% 
  select(-label)
mtb_standard_scaled <- scale(mtb_numeric)

mtb_numeric <- mtb_no_null %>% 
  select(-label)

mtb_numeric <- mtb_no_null %>% 
  select(-label)

clusters <- 1:10
dists <- c()
for (c in 1:10){
  km <- kmeans(mtb_standard_scaled, centers=c, iter.max=1000)
  dists <- c(dists, km$tot.withinss)
}

# jpeg('../Images/Kmeans.jpg')
# plot(clusters, dists, type='l', xlab='Clusters', ylab='Total Sum of Squared Euclidean Distances')

# Plot our results
tibble(clusters = clusters,
       dists = dists) %>% 
  ggplot(aes(x = clusters, y = dists)) + 
  geom_point(size = 3, alpha = .9, color = 'slateblue') +
  geom_line(size = 2, alpha = 1, color = 'slateblue1') + 
  labs(title = "K-Means Clustering of MTB Data",
       subtitle = 'Method uses `tot.withinss` parameter to measure distances.',
       x = 'Clusters',
       y = 'Total Sum of Squared Euclidean Distances')

# Let's see where these clusters would end up on the 2D PCA plot


mtb_pca_scaled <- prcomp(mtb_standard_scaled,
                  center = F,
                  scale. = F)

pca_2_scaled <- as.matrix(mtb_standard_scaled) %*% as.matrix(mtb_pca_scaled$rotation[,1:2])

pca_km_scaled <- kmeans(pca_2_scaled, centers=3, iter.max=1000)

#Something's not working here, definitely hitting a local min or something
colorgroups <- function(g){
  if (g == 'tr' || g == 'Trail'){
    return('blue')
  }
  else if (g == 'xc' || g == 'Cross Country'){
    return('pink')
  }
  else if (g == 'dc' || g == 'Downcountry'){
    return('darkgoldenrod3')
  }
  else if (g == 'am' || g == 'All Mountain'){
    return('red')
  }
  else if (g == 'en' || g == 'Enduro'){
    return('green')
  }
}
catNames <- c('Cross Country', 'Downcountry', 'Trail', 'All Mountain', 'Enduro')
cols <- unlist(lapply(new_mtb_data$label, colorgroups))

unlist(lapply(catNames, colorgroups))

# jpeg('../Images/PCA_clusters.jpg')
plot(pca_2_scaled, col=cols)
points(pca_km_scaled$centers, col = 'slateblue4', pch = 'x', cex = 1.5)
# text(pca_2_scaled[,1], pca_2_scaled[,2], rownames(pca_2_scaled))
legend("bottomleft", legend= c(catNames, 'Cluster Center'), col=c(unlist(lapply(catNames, colorgroups)), 'black'), pch=c(rep('o', 5), 'X'))

#TODO let's look at this bottom cluster - both Niner bikes
#Niner has low reach numbers on its bikes - could be because we used the Medium for these!
#Based on PCA mapping, the blur tr, expic, Exie, Ripley, and Element all have less chainstay length, and less pistons?? wow, should we exclude piston count?? with more 2 piston bikes getting added, it evens out the average, so these aren't showing up as much anymore

p_load(ClusterR)

# Build our GMM model
mtb_gmm <- GMM(mtb_standard_scaled,
               dist_mode = 'eucl_dist', # Distance metric to use during seeding of initial means clustering
               seed_mode = 'random_subset', # How initial means are seeded prior to EM alg
               km_iter = 10, # Num of iterations of K-Means alg
               em_iter = 10, # Num of iterations of EM alg
               verbose = T
               )

mtb_gmm_pred <- predict(mtb_gmm, mtb_standard_scaled)

opt_gmm <- Optimal_Clusters_GMM(mtb_standard_scaled, 
                               max_clusters = 20, 
                               criterion = "BIC", 
                               dist_mode = "eucl_dist", 
                               seed_mode = "random_subset",
                               km_iter = 10, 
                               em_iter = 10, 
                               var_floor = 1e-10, 
                               plot_data = T)

p_load(mclust)

mtb_gmm2 <- Mclust(mtb_standard_scaled)

#or specify number of clusters
# mb3 = Mclust(iris[,-5], 3)

# optimal selected model
# mtb_gmm2$modelName

# optimal number of cluster
# mtb_gmm2$G

# probality for an observation to be in a given cluster
# head(mtb_gmm2)

# get probabilities, means, variances
summary(mtb_gmm2, parameters = TRUE)

plot(mtb_gmm2, 'classification')

p_load(e1071,
       caret)

#convert all mountain category to enduro, dc -> Xc?
remap <- function(x, num){
  if (x=='am' || x=='en'){
    if (num){
      return(4)
    }
    else{
      return('Enduro')
    }
  }
  else if(x=='xc'){
    if(num){
      return(1)
    }
    else{
      return('Cross Country')
    }
  }
  else if(x=='dc'){
    if(num){
      return(2)
    }
    else{
      return('Downcountry')
    }
  }
  else if(x=='tr'){
    if(num){
      return(3)
    }
    else{
      return('Trail')
    }
  }
}
labels <- as.factor(unlist(lapply(new_mtb_data$label, remap, F)))
n <- length(labels)
test_idx <- sort(sample(1:n, round(n/5)))

Xtest <- mtb_standard_scaled[test_idx, ]
Xtrain <-  mtb_standard_scaled[-test_idx, ]

trainSVM <- function(x, y, idx){
  xtest <- x[idx,]
  xtrain <- x[-idx,]
  ytest <- y[idx]
  ytrain <- y[-idx]
  
  clf <- svm(x=xtrain, y=ytrain)
  preds <- predict(clf, xtest)
  
  acc <- 0
  cm <- table(ytest, preds)
  for (i in 1:length(unique(labels))){
    acc <- acc + cm[i,i]
  }
  return(acc/sum(cm))
}

folds <- createFolds(labels, k=10)

accs <- c()

for (fold in folds){
  acc <- trainSVM(mtb_standard_scaled, labels, fold)
  accs <- c(accs, acc)
}

mean(accs)

# Roughly 60% accuracy when treating down country as separate category
# But - roughly 77% accuracy when treating down country as XC, only 65% accuracy when treating downcountry as trail, suggests that downcountry bikes are more akin to XC than they are trail

folds <- createFolds(labels, k=10)

accs_2pc <- c()

for (fold in folds){
  acc <- trainSVM(pca_2_scaled, labels, fold)
  accs_2pc <- c(accs, acc)
}

mean(accs_2pc)

dat <- data.frame(x=cbind(pca_2_scaled[,2], pca_2_scaled[,1]), y=labels)
pcsvm <- svm(y~., data=dat)

# jpeg('roughSVM.jpg')
plot(pcsvm, dat)

# circle -> correctly predicted
# X -> incorrectly predicted
# black -> true XC
# Red -> true DC
# Blue -> true TR
# Green -> true EN